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Abstract. Concentrations of matter, such as galaxies and galactic clusters, 
originated as very small density fluctuations in the early universe. The existence 
of galaxy clusters and super-clusters suggests that a natural scale for the matter 
distribution may not exist. A point of controversy is whether the distribution is 
fractal and, if so, over what range of scales. One-dimensional models demonstrate 
that the important dynamics for cluster formation occur in the position- velocity plane. 
Here the development of scaling behavior and multifractal geometry is investigated for 
a family of one-dimensional models for three different, scale-free, initial conditions. 
The methodology employed includes: 1) The derivation of explicit solutions for the 
gravitational potential and field for a one-dimensional system with periodic boundary 
conditions (Ewald sums for one dimension); 2) The development of a procedure for 
obtaining scale-free initial conditions for the growing mode in phase space for an 
arbitrary power-law index; 3) The evaluation of power spectra, correlation functions, 
and generalized fractal dimensions at different stages of the system evolution. It is 
shown that a simple analytic representation of the power spectra captures the main 
features of the evolution, including the correct time dependence of the crossover from 
the linear to nonlinear regime and the transition from regular to fractal geometry. A 
possible physical mechanism for understanding the self-similar evolution is introduced. 
It is shown that hierarchical cluster formation depends both on the model and the 
initial power spectrum. Under special circumstances a simple relation between the 
power spectrum, correlation function, and correlation dimension in the highly nonlinear 
regime is confirmed. 



1. Introduction 

The observation that visible matter in the universe shows structure on a huge range 
of scales, from galaxies, to clusters, to super-clusters, to voids pQ, has led Mandelbrot, 
Pietronero and others to conjecture that the structure of the universe may be fractal 
[2]. Support for this controversial conjecture is provided by the fact that the inter- 
galactic, two-body correlation function decays as a power-law [3]. If the geometry of 
the universe was truly fractal on all scales, then the basic principle of cosmology, that 
the universe is homogeneous and isotropic on large scales, would be violated. While 
current observations support the existence of an upper bound for the size of the largest 
structures, the issue is not completely closed [H [3]. Typically, in nature, observable 
fractal behavior is usually restricted to a finite scale range [6]. Regardless of whether 
such a bound exists, it is still important to understand how the appearance of fractal 
geometry develops from the underlying dynamics. 

The remarkable uniformity of the temperature distribution of the cosmic 
background radiation CMB is a consequence of the nearly uniform energy distribution at 
recombination [7]. Over time, the small (order 10~^) density fluctuations developed into 
the structures we see today. Measurements and theory suggest that, following inflation, 
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the density fluctuations were Gaussian random variables with an approximate power- 
law spectral density. They were subsequently modulated by baryonic acoustic oscillation 
before recombination [3] . The commencement of structure formation in the dark matter 
component actually preceded recombination. These circumstances provide a natural 
scenario for the initial conditions of any dynamical study. With current technology 
astronomers can observe bright objects in the distant past. However, they cannot 
directly observe dynamical processes that take giga-years to unfold. As a consequence, 
computer simulation plays an especially prominent role in astrophysics, and has been 
employed to investigate complex behavior ranging from galaxy collisions and mergers 
to structure formationjS]. Although progress has been substantial, the ability of three- 
dimensional simulations with on the order of 10^ particles to resolve fractal scaling laws 
is still hampered by limited resolution and approximations in the underlying dynamics. 
In contrast, simulations with one-dimensional models can incorporate up to 10^ particles 
without compromising the two-body gravitational interaction. Thus, although they are 
less realistic "toy" models, they have the potential to yield insights into clustering with 
current computers. In addition, since theory suggests that the first collapsed objects 
were highly flattened "pancakes", they may yet provide contact with the real universe 

EE]- 

Rouet and Feix were the first to recognize the potential for using one-dimensional 
models to investigate clustering in a matter-dominated cosmological setting [TOl [TT] . 
Since there is no curvature in one dimension, general relativity does not provide a 
unique path for obtaining the correct laws of motion. Nonetheless, they showed how 
the transformation to co-moving coordinates could be accomplished in a completely 
self-consistent formulation which we refer to as the RF (Rouet-Feix) model. Starting 
with a spatially uniform initial distribution, they demonstrated that, as time evolved, 
hierarchical clustering occurs. Thus the one-dimensional evolution is similar to what 
is believed to have occurred in the universe following recombination. Rouet, Feix and 
Jamin computed the box counting dimension for the distribution of matter in the one- 
dimensional configuration and /i (position, velocity) space. They obtained a result of 
about 0.6 for the configuration space, strongly suggesting a fractal geometry. Since 
their seminal work, another similar model was proposed by Fanelli and Aurell, referred 
to simply as the quintic or Q model [12]. As we will see, the Q model sacrifices internal 
consistency in order to maintain the correct coefficient of the average cosmological 
density. 

A different approach employing a one-dimensional model to explore matter- 
dominated cosmology was taken by Gouda and co-workers. Yano and Gouda employed 
the Zeldovich approximation [9] to investigate the evolution of the one-dimensional 
system [13]. When there is no cutoff in the initial spectra they showed that the evolution 
of the power spectra is self-similar, and established that three different scaling regimes 
occur, each with its own characteristic spectral index. Later Yano et al [H] studied 
the development of a single wave in the phase plane to investigate the effect of caustics 
on the evolution of the power spectra. In addition, Tatekawa and Meida used the 
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Zeldovich approximation to study the self-similar evolution of a one-dimensional system 
with initial conditions selected from a Cantor set [15]. 

More recently Miller and Rouet extended the original work of Rouet and Feix 
to include an investigation of multi-fractal properties [HI [17]. They considered the 
fact that the clusters are actually forming in /i space. An analysis of the generalized 
fractal dimensions was applied to both configuration and /i space and it was shown that 
multifractal geometry develops in each. Miller and Rouet carried out a multifractal 
analysis of both the RF and Q models, as well as a Hamiltonian version. Valageas has 
investigated the Hamiltonian version with identical boundary conditions, and has shown 
the existence of a sequence of stationary and equilibrium states [IHl 112]. The use of one- 
dimensional models has recently been extended by Sutter and Ricker to include a dark 
energy field in addition to dark matter [20] • They investigated how the introduction 
of dark energy influences the formation of Zeldovich "pancakes", presumably the first 
large-scale structures to emerge in the cosmos. In a recent series of articles Gabrielli et al 
have studied the one- dimensional system in the infinite particle limit [211 [22l [23l [2H EH] . 
A central focus of the series is the statistical properties of the force distribution obtained 
from sampling static distributions of particle positions from infinite perturbed lattices. 
In [23] a number of the issues mentioned here are discussed within this framework. 

Theoretical cosmology suggests that during the period that density fluctuations 
remain small, evolution is linear and is dominated by a "growing mode" and the spectral 
distribution of the fluctuations remains a power- law [1]. However, in the simulations 
of Miller and Rouet described above, initial conditions were created by independently 
sampling the velocity of each particle from either a uniform or Gaussian distribution. In 
the work presented here, the initial conditions follow the cosmological picture for cold 
dark matter. We have performed simulations and a multi-fractal analysis of each model 
for a variety of initial, power-law, spectral densities. Here we describe the circumstances 
under which the evolution follows the accepted pattern of self-similar, hierarchical 
clustering. We show how the evolution of the spectral density and correlation function 
depend on both the model and the initial power-law index. We find that, for sufficiently 
large samples, the analytic relation between correlation dimension, correlation function, 
and spectral density is obeyed. 

In the following we first explain the mathematical formulation of the different 
models in section [2j In particular, we show how to construct the exact gravitational 
potential and field for a one- dimensional, periodic system, i.e., Ewald sums for one 
dimension. In Section |3] we describe how simulations are performed and show how to 
create scale-free initial conditions for the growing mode in phase space for an arbitrary 
initial scaling index. We also present typical results for the Q model in the highly 
nonlinear region, including the calculation of power spectra and the two-body correlation 
function. In addition a simple analytic representation of the power spectra is introduced 
that captures the main features of the evolution, including the correct time dependence 
of the crossover from the linear to nonlinear regime and the transition from regular 
to fractal geometry. In Section |4] we introduce techniques for measuring generalized 
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fractal dimensions and apply them to the Q model. In addition we give a comparison of 
results for related models and different initial power spectra and demonstrate that the 
power spectra, correlation function, and correlation dimension are interrelated in the 
nonlinear regime. To provide a physical basis for the observed behavior we introduce the 
mechanism of hierarchical virialization. Finally, in the last section, we present general 
conclusions that can be drawn from this investigation and provide a discussion of some 
open questions that still need to be addressed. 

2. One-Dimensional Systems 

Although dynamics in the universe is governed by the general theory of relativity, in 
a sufficiently small sample Newtonian dynamics provides an adequate working model 
In order to study structure formation, cosmologists focus on a segment of the 
universe Q that is small enough that a Newtonian description is adequate, but much 
larger than the two-body correlation length p]. They represent the location of the mass 
points of N-body simulations in comoving coordinates that follow the Hubble flow, so 
the average density remains constant. They ascribe a cubical shape to the sample and 
assert periodic boundary conditions to obtain a closed dynamical system that mimics 
a segment of the real universe [27]. This conserves mass and insures the continuity 
and smoothness of the gravitational potential and field at the boundaries. A potential 
problem is that, on the average, the universe is isotropic but the symmetry imposed by 
the cubical boundary conditions is not. For a large enough cube, this problem can be 
avoided for a finite time. In tree models Ewald sums are typically employed to compute 
the gravitational force from the infinite number of system "copies" [28] while in grid 
models the periodic boundary conditions are built into the potential [2S]- 

A one-dimensional gravitational system corresponds to a set of parallel mass sheets 
moving in the direction perpendicular to their surface. In the seminal work of Rouet 
and Feix, they assumed a universe with spherical symmetry about a point. The 
elements of their universe were then concentric, spherical, irrotational mass shells. 
Far from the symmetry center, the radius of curvature is large and, locally, the shells 
are approximately planar and parallel. By choosing such a segment, they obtained a 
stratified system of comoving, planar, mass sheets for their model. However, if the 
object is to construct a consistent, one-dimensional, gravitational model that embraces 
an expanding background, it is not necessary to assume any special symmetry. One can 
just start with one dimension and assert a uniform expansion factor that also applies 
to the parallel dimensions of the mass planes. If gravity is the only force acting, then 
internal consistency forces the t^/^ time dependence associated with the expansion factor 
of the Einstein de-Sitter, or matter-dominated, universe. However, the coefficient of 
the time-dependent mean density differs from the regular, three-dimensional model. 
To avoid this mild conundrum Fanelli and Aurell took a different approach. They 
first transformed to comoving coordinates in homogenous and isotropic three-space and 
then inserted a system of planes. Because the negative background density induced 
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by the transformation to comoving coordinates doesn't cancel with the average density 
associated with the mass sheets, an additional correction has to be included [17]. Here 
we will follow this approach. 

Consider such a bounded region Q. We are interested in the evolution of density 
fluctuations following the time of recombination, so that electromagnetic forces can be 
ignored and Newtonian dynamics provides an adequate representation of the motion in 
a finite region Then, in a (3+l)-dimensional universe, the Newtonian equations 
governing a mass point are simply 

c^i" (iv ^ , , , . 

-=v. ^=E,M (1) 

where, here, Fig{r,t) is the gravitational field. To follow the motion in a frame of 
reference where the average density remains constant, i.e. the comoving frame, we 
introduce the scale factor A{t) for a matter-dominated universe p] and transform to a 
new space coordinate which scales the distance according to A{t). Writing r = y4(t)x 
we obtain 

(i^x 2 c/A (ix 1 (PA 1 

1^ + Ali^ + Al^"" =^E,(x,t) (2) 

where, in the above, we have taken advantage of the inverse square dependence of the 
gravitational field to write Eg(x, t) = ^Eg(r,t) where the functional dependence is 
preserved. In a matter-dominated (Einstein-de Sitter) universe we find that 

2 

A{t)=l^-y, p,it) = {enGt^y (3) 

where to is some arbitrary initial time corresponding, say, to the epoch of recombination, 
G is the universal gravitational constant, and pb{t) is the average, uniform, density 
frequently referred to as the background density. The justification for equation ^ 
comes from the Robert son- Walker metric and the Friedman equation [I] . However, these 
results can also be obtained directly from equation ^ by noting that if the density is 
uniform so that all matter is moving with the Hubble flow, the first two terms in equation 
(|2| vanish whereas the third term (times A) must be equated to the gravitational field 
resulting from the uniformly distributed mass contained within a sphere of radius A |x| 
which is simply the right-hand side of equation (|2|. Then the third term of equation ^ 
is the contribution arising by subtracting the field due to the background density from 
the uniform sphere [T]. Noting that A^ph{t) = pb{to) forces the result. Alternatively, 
also for the case of uniform density, taking the divergence of each side of equation ^ 
and asserting the Poisson equation forces the same result. Thus the Friedman scaling 
is consistent with the coupling of a uniform Hubble flow with Newtonian dynamics [1]. 

In standard three-dimensional cosmological simulations it is common, but not 
ubiquitous, to employ conformal time, i.e to use A(t) as the measure of progress [2^ . 
Here, for computational purposes, we will see that it is useful to obtain autonomous 
equations of motion with coefficients that do not depend explicitly on the time. This 
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can be effectively accomplished flU[ [TT] by transforming the time coordinate according 
to 

dt = B{t)dT, B{t) = - (4) 

^0 



yielding the autonomous equations 
rf^x 1 dx 2 



x = E,(x). (5) 



rfr2 3to dr 9t§ 

These can be further simplified by choosing the inverse Jeans' frequency Tj for the unit 
of time [30] 



Tj = uj-/ = {A'KGpr'/' = ^^-t, (6) 

yielding 

c?^x 1 (ix 1 , . 

where we used the fact that 3tQ/2 = 1 in the adopted time units and r is now expressed 
in the new units. Thus equation ^ corresponds to a dissipative dynamical system in the 
comoving frame with friction constant 1/a/6 and with forces arising from fluctuations 
in the local density with respect to a uniform, three-dimensional, isotropic, neutralizing 
background. 

We now imagine that the actual source of the density fluctuations is a system 
of parallel mass sheets with a neutralizing background, similar to a single- component 
plasma [311 E2] • For the special case of the stratified mass distribution induced by the 
one-dimensional system, the local particle density at time t is given by 

p{x, ^) = X! — Xj) (8) 

where mj{t) is the mass per unit area of the j*'* sheet and, from symmetry, the 
gravitational field only has a component in the x direction. To avoid confusion, we 
will just refer to this as the mass. From Gauss' law we know that the field due to a 
single mass sheet is constant. Then, for an isolated system with a finite number of 
particles and no background, the field experienced by one of the sheets is proportional 
to the difference between the mass on each side [33] • To obtain equations of motion, if 
we just take the component of equation ([T]) in the x direction we see that a problem 
arises. The third term corresponds to the field from a sphere of constant density, but 
we now have a slab geometry so this term needs to be multiplied by a factor of three to 
maintain mass neutrality [T7] . 

Following the three dimensional example, we will also assume periodic boundary 
conditions to approximate the actual behavior in a finite slab of our model one- 
dimensional universe for a finite time. The size of the system will be determined by the 
length of elapsed time we require to evolve the system before the presence of boundaries 
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play a significant role. To accomplish this we need to determine the gravitational field 
induced by a single particle in the periodic system. Assume that the width, and hence 
the spatial period, of our system in the comoving frame is 2L, so we may choose the 
position xi of our particle of mass mi in [—L, L) with the points at L and —L identified. 
For periodic boundary conditions the potential can only be defined for a mass-neutral 
system [291 EZl EH] where the source of the potential and field is the difference between 
the local density and its average over one period. Since our particle carries with it 
the negative background density — mi/2L, the potential it induces, 0i(a;), satisfies the 
following Poisson equation 

= 4vrmiG[5(x - x,) - ^] (9) 
with the general solution 

01 (x) = 27rmiG'[|x — Xi\ -(x — Xif + h{x — Xi) + c] (10) 

where h and c are constants. 

The periodicity requirement, 0i(L) = 0i(— L), forces 6 = 0. In addition, to 
guarantee that as L — )■ oo the potential approaches that of the isolated system, we 
choose the arbitrary additive constant c = 0, yielding for the final form of the potential 

01 (x) = 27rmiG[|x — xi| — — xi)^]. (11) 

We immediately obtain the gravitational field Ei induced by a single particle in the 
periodic system: 



Eiix) = = 27rmiG 

ox 



— (x — Xi) + 0(xi — x) — 0(x — Xi) 



(12) 



where 0(x) is the usual step function and, for consistency, we assign 0(0) = |. Note that 
we could have obtained the result from symmetry as well. Periodic boundary conditions 
in one dimension restrict the motion to the one-torus or circle. Since there is no preferred 
direction on the circle, the spatially averaged field must vanish. Thus symmetry alone 
guarantees that 6 = 0. The results for the potential and field, equations ( lTp2 ) can also 
be obtained by employing a screening function to directly sum over the periodic system 
images, or replicas [31], and are effectively Ewald sums for one dimension. Unlike their 
three-dimensional analogs [251 , they have a simple closed-form expression. 

Let us assume that the primitive cell of our periodic system contains 2N particles 
(sheets) confined within a slab with width 2L, i.e. —L < x < L. For the special case of 
equal masses rnj{t) = m{t), by direct summation the correct form of the gravitational 
field occurring at the location of particle i in comoving coordinates is then 

2N 

Eg{xi) = 27im{to)G[—{xi - xj + Nr^^ - NL,i] (13) 

since we already implicitly accounted for the fact that mj{t) = m(to)/^^ in equation 
(|5|. In equation (13) Nr^i {NL,i) is the number of sheets on the right (left) of particle 
i and Xc(t) is the system center of mass in [— L, L). The field is further simplified by 
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establishing the connection between m{to) and the background density at the initial 
time, Pb(to). Then 

pM = {GnGtiy' = (^j^ m(to), (14) 
and we may express the field by 

Eg{x,) = ^ (^^) [^(x, - + Nn, - Nl,]. (15) 



(T^ [2{xi - Xc{t)) + Nn,^ - NL,i]. (16) 



The equations of motion for a particle in the system now read 
dr'^ ^ y/Q dr 

where we have taken account of the adopted unit of time (see equation ([6])) and chosen 
the mean inter-particle spacing L/N for the unit of length. 

The description is completed by noting that, since the system satisfies periodic 
boundary conditions on the interval [—L, L), when a particle leaves the primitive cell 
defined by — L < x < L on the right (left), it re-enters at the left (right) hand boundary 
with the identical velocity. Note that from equation (13) there is no change in the field 
experienced by other particles during such a boundary crossing. This is a necessary 
condition for the one-torus geometry, since there is nothing special about this point. 
Other consequences of the torus geometry are discussed in [31]. In our simulations 
we made the further assumption of symmetry about the origin: For each particle with 
< X < L with velocity v, there is a twin located at —x with velocity —v. With this 
stipulation, Xc = in equation ( 16 ) and the periodic boundary conditions are equivalent 
to an A^-particle system with refiecting boundary conditions at x = 0, L. For the times 
of interest, we demonstrated in [M] that the behavior of the system is nearly identical for 
each type of boundary. Note that, without the symmetry requirement, if one arbitrarily 
sets Xc = the torus geometry is violated. 

As we mentioned earlier, by embedding the stratified system into a region of three- 
dimensional Euclidean space we are changing the local symmetry, so we have to adjust 
the local background density. In concept it is similar to the well-known Zeldovich 
approximation that is used to investigate the formation of the first matter concentrations 
that are believed to have a "pancake" geometry [1] . In their investigation of the structure 
of Zeldovich pancakes Aurell et al have shown that the models are closely related [36] . 
In [23] it is argued that they are equivalent. Notice that this hybrid model mixes two 
different symmetries, the isotropic 3+1 dimensional system that cosmologists model with 
a periodic cube, and the planar geometry that is represented by a slab. Since the mass 
concentration in the idealized planes stretches to infinity in two directions it cannot 
represent the physical universe. The slab geometry can only represent cosmological 
evolution for a short time before the infiuence of the other dimensions (terms in the 
strain tensor) become important [1]. In contrast, the RF model is obtained from the 
reverse sequence where one first restricts the geometry to 1+1 dimensions and then 
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introduces the transformation to the comoving frame. In this approach the derivation 
is more straightforward; it is not necessary to make the adjustment in the coefficient 
of Xi as we did here to obtain the correct background contribution [ini [HI HE] . This is 
quickly seen by noting that the divergence of x is three times greater than the divergence 
of XX which one would obtain by directly starting with the one-dimensional model. 
As a consequence, in the RF model, the coefficient of the first derivative term (the 



friction constant) in equation (16) is l/\/2 instead of In every other respect 

the models are identical. This simply illustrates that, since there is no curvature in 
a (l + l)-dimensional universe, there is a degree of arbitrariness in choosing the final 
model. It cannot be obtained solely from general relativity. For a discussion of this 
point, see Mann et al [37j. A Hamiltonian version can also be considered by setting the 
friction constant equal to zero. In their earlier work, using the linearized Vlasov-Poisson 
equations, Rouet and Feix carried out a stability analysis of the model without friction. 
They determined that the system followed the expected behavior, i.e., when the system 
size is greater than the Jeans' length, instability occurs and clustering becomes possible 
flU\ [TT] . Of course, when the friction term is not present, both the Q and RF models 
are identical. Then, with the assumption that the friction term will not have a large 
infiuence on short-time linear stability, the analysis of the Hamiltonian version applies 
equally to both versions. 

3. Simulations and Initial Conditions 

An appeal of these one-dimensional gravitational systems is their ease of simulation. 
In each of the one-dimensional Newtonian models considered here (Q, RF and 
Hamiltonian), it is possible to analytically integrate the motion of the individual 
particles between crossings. Then the temporal evolution of the system can be obtained 
by following the successive crossings of the individual, adjacent, particle trajectories. In 
particular, for the Q model, if we let yi = Xj+i — Xi, where we have ordered the particle 
labels from left to right, then we find that the differential equation for each yi is the 
same, namely 

d^Vi , ^ dyi 



The general solution of the homogeneous version of equation (17) is a sum of 
exponentials. By including the particular solution of the inhomogeneous equation 
(simply a constant) we obtain a fifth order algebraic equation in m = exp(r/\/6) for 
the successive crossings, defined by yjir) = ; hence the name Q, or quintic, model. 
These can be determined numerically in terms of the initial conditions by analytically 
bounding the roots and employing a numerical root-finding method. Note that for the 
RF model a cubic equation is obtained so the crossing times can be found analytically 
[Tot [TTl [16] . A sophisticated, event-driven algorithm was designed to execute the 
simulations. Using the Newton-Raphson method, the algorithm computes all possible 
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crossings of adjacent pairs of particles with double-precision accuracy. Since, between 
crossings, the general solution of the dynamical equations is known analytically, once a 
crossing time is estabhshed, the position and velocity of each particle can be determined 
with the accuracy of the computer. Two important features of the algorithm are that it 
only updates the positions and velocities of a pair of particles when they actually cross, 
and it maintains the correct ordering of each particle's position on the line. In contrast 
with typical N-body simulations, it is not necessary to introduce a discreteness length or 
time step. Using the algorithm we are able to carry out runs for significant cosmological 
times with large numbers of particles. In particular, it is possible to simulate a system 
until all of the mass is concentrated in just a few clusters. Depending on the initial 
conditions, this typically occurs on the order of 15-20 dimensionless time units into the 
simulation. Since, at this stage, the influence of the boundary conditions can no longer 
be ignored, there is no advantage in continuing the runs any further. 

In cosmology, a consequence of inflation is that primordial density fluctuations are 
independent Gaussian random variables with a scale-free, power-law, spectral power 
distribution P{k) ~ A;" [1]. Analysis of the WMAP measurements of the angular 
distribution of the CMB [7| yields a value for rig very close to unity, the Harrison- 
Zeldovich value. In this study we have investigated the dynamics of all three models 
with three different initial power spectra with indices, n = 0, 1, 2. To construct the 
initial conditions we assume that the system is nearly uniform. If we use an ordered 
labeling of the particles, then the equilibrium positions are simply Xjo = ^(j — |) for 
j = 1, ...,N. Initially the system is in the linear regime and no crossings have occurred 
so the particles are assumed to be very close to their equilibrium points. Then, from 
equation (13) it is straightforward to show that the gravitational field on particle j 
with position xj arising from the initial density fluctuation is, in our units, simply 
SE{xj) Represent the local density fluctuation 5p{x) as a Fourier series 

with coefficients 5pk- Then, from the Poisson equation, 

5E{x,) = -ATxGY,^5pu{e^^^ - 1) 

k ^ 

where we note that, at the origin, 5E vanishes and there is no contribution from = 0. 
Taking the Fourier coefficients as Gaussian random variables with a power-law power 
spectrum \6pk^ of index n, we have 5pk ~ /^"/^g^^fc where the random phases dk are 
uniformly distributed on [0, 27r] and 9-k=—9k to insure 6p is real. Combining the above 
we obtain 

Zj = Xj — Xjo = C ^ k"'^'^^^[sin{kxjo + 6k) — sin{6k)], k = vtt/N, v = 1 

k>0 

for the initial displacements of the particles from equilibrium. The sum is cut off 
when the wavelength is smaller than the mean inter-particle spacing. The constant of 
proportionality C was chosen just small enough to prevent any crossings of the ordered 
positions. In the linear regime, from equation (16), the evolution of each particle is 
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governed by 

d'^Zj 1 dz. 



+ = ^E{x,) = z, (18) 




cir2 ^ dr 
with general solution 

Zj{T) = aexp ^-^r^ + bexp 

Following standard practice, we include only the contribution from the growing mode. 
This fixes the initial velocity of particle j as -^zj and completes the assignment of the 
initial condition of the N-body system. For an alternative approach for indices n = 0, 1, 2 
see [H]. 

In figure [l] we present a visualization of a typical run with 2^^ particles. The system 
consists of the Q model with an initial spectral index of n = 2. Initially the velocity 
spread is small, within (-1, 1) in the dimensionless units employed here, and the system 
contains on the order of 10^ Jeans' lengths. Of course, since the initial state does not 
represent thermal equilibrium, the Jeans' length lacks predictive certainty, but strongly 
suggests that instability will ensue. Following the original work of Rouet et al fiU\ [TT] . 
we present a sequence of snapshots at integer values of the dimensionless scaled time 
through T = 18. In the left column we present a histogram of the particle positions 
at increasing time frames, while on the right we display the corresponding particle 
locations in /x space. It is clear from the panels that hierarchical clustering is occurring, 
i.e., small clusters are joining together to form larger ones, so the clustering mechanism 
is "bottom-up" [1]. The first clusters seem to appear at about r = 8 and there are many, 
while by r = 18 there are only on the order of a dozen clusters. In larger simulations 
in the space we observe that, between the clusters, matter is distributed along linear 
paths [17]. As time progresses the size of the linear segments bridging the clusters 
increases. The behavior of these under-dense regions is governed by the stretching in 
H space predicted analytically by Vlasov theory [161 IE!- Qualitatively similar histories 
are obtained for the RF model but the clustering occurs more slowly. However, there 
are some subtle differences. In figures [2] we zoom in and show magnified inserts from 
the mass distribution in /i space at time r = 18. The hierarchical structure observed in 
these models suggests the existence of fractal geometry, but careful analysis is required 
to determine if this is correct. 

In figure [s] we show the power spectra of the initial density fiuctuations |5pfc|^. 
Although the data are noisy, by averaging over a group of neighboring points it is 
easy to extract the slope with good accuracy. Note that the graph conforms with 
the construction of the initial state described above. This confirms the validity of the 
initialization procedure. Note also that on small scales {k > 10) the slope fiattens to 
showing that, for wavelengths less than the inter-particle spacing, the power spectra is 
simply a white noise as anticipated. 

Historically, power-law behavior in the tail of the density-density correlation 
function has been taken as the most important signature of self-similar behavior of the 
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Figure 1: Evolution in configuration and fi space for tlie quintic model with 2^^ particles 
from r = to r = 18. The initial distribution is such that the density power spectrum 
has a power-law of index n = 2 (see Figure Isj). 



(b) (c) 

Figure 2: Consecutive expansions (zooms) in the fi space panels. Panel (a) shows the 
region selected by a rectangle in the /i space panels at r = 18 of figure [l} Panel (b) is 
the region selected by a rectangle in panel (a) and so on. These representations have 
the appearance of a random fractal which suggests self-similarity. 



distribution of galaxy positions [3l [38] . Even at r = 18 we can see from the snapshots 
that the clusters are well developed, but still much smaller than the system size. Then 
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Figure 3: The Power spectrum of the initial density fluctuations for the quintic model. 
Particles have been distributed so that the power spectrum presents a power-law of 
index n = 2. 



it is safe to assume that boundary conditions do not yet play a significant role. In figure 
4] we provide a log-log plot of the correlation function ^(r) at r = 13 defined by 

^{r) = {6p{x + r)6p{x)) 

r+A 



1 

A 



J2 Sir' - (Xi-Xj)) - 



N{N -1) 
2Z2 



dr' 



(19) 



where particles i and j are such that L/4 < Xj < 3L/4 to avoid boundary effects, and A 
is the bin size. We observe three distinct regions in figure |4j First there is a relatively 
fiat region at small scales. It is followed by a scaling region from about 0.02 < r < 2.0, 
a range of about 2 decades, where ^(r) oc r"''', which finally deteriorates into noise. 
This is similar to cosmological studies which also exhibit the intra-cluster structure at 
small scales, and the inter-cluster scaling seen here. In addition, in the real cosmological 
setting, there is a peak due to baryonic acoustic oscillations [3]. The appearance of noise 
at larger scales occurs because the correlation function is decreasing while the noise is 
relatively constant. In a log- log plot, as r is increased, at some point the noise becomes 
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as large as the mean so large fluctuations appear in the data. Expressed differently, in a 
log-log plot of the correlation function, the apparent noise is effectively magnified with 
increasing r until it dominates the mean. A straightforward linearization shows that 
the effective amplitude of the noise in InC grows as exp(7lnr) while it is still relatively 
small. In principle it should be possible to reduce the noise by taking an ensemble 
average of ^(r) over many runs, and then computing the slope of ln(^(r)). However, 
since we already employ 2^^ particles, and the effective noise is growing exponentially 
with increasing r, this may only slightly extend the apparently noise-free region. By 
computing the slope 7 of the log-log plot in figure |4] we are able to obtain an estimate 
of the correlation dimension D2 = 1 — '~f ioi a one-dimensional system [52] • We find 
that 7 = 0.373 for a scahng region of about two decades in /. This suggests that the 
correlation dimension is approximately 0.63, which is in agreement within the standard 
numerical error with the multifractal analysis described below in some detail (see table 




Figure 4: The correlation function at r = 13 for the quintic model. It exhibits a 
scaling region in r from about 0.02 to about 2, a range of about 2 decades with a scaling 
exponent 7 = .373. 



Cosmology in One Dimension: Fractal Geometry, Power Spectra and Correlation 16 



Table 1: Exponent values for the three different models. For pure power law behavior 
D2 = 1 — 7 = —n. For each system, = 65535 and, initially, n = 2. 



Model 


Time T 


7 


n 


D2 


Q 


13 


.373 


-.672 


.623 




16 


.349 


-.637 


.628 


RF 


13 


.552 


-.490 


.478 




16 


.547 


-.496 


.500 


H 


13 


.173 


-1.048 


.862 




16 


.163 


-.880 


.853 



We have followed the evolution of the power spectra of the density fluctuations 
in these systems. In figure |5] we display log-log plots of P{k) for a sequence of times. 
We see from the plot that two scaling regimes are present. For small k, i.e. on large 
scales, the system appears to retain the same dependence on k as in the initial state 
and the linear regime. In contrast, using the same averaging method as earlier, we 
find that a second power law dependence develops at large k. At r = 13 , for k less 
than .02 to more than 50, i.e. over a range greater than three decades, this power-law 
behavior is clearly exhibited. As time progresses the division between the two regimes 
moves to the left so the shape is preserved and the process is self-similar. Between the 
two scaling regions is a transition region centered at, say, k = kdr) which we define 
below. By r = 16 there is no longer any trace of the linear regime. The behavior 
for k < kc is a. remnant of the inter-cluster regions where the particles have not yet 
crossed. From Eq. ( 18 ) and the following discussion we see that the initial amplitude of 
deviations from the particle equilibrium positions in the growing mode is increasing as 
exp(2r/-\/6) resulting in the corresponding power spectra in the linear region to increase 
as exp(4r/-\/6). Analysis of the plots in Fig. (jsj) confirm this nicely. This amplification 
corresponds to the "stretching" in the low density regions of fi space predicted by Vlasov 
theory [121 [IZj. It is well known, and discussed from a fluid picture elsewhere in the 
literature [HllQlES]. 

Notice that in Fig. ([s]), as a rough approximation, the power spectra in the 
nonlinear region has the form Pni{k) = Bk^ where B is approximately constant and 
n ~ —.65. Therefore, in both the linear and nonlinear regimes, the power spectra 
are well characterized. We can use this rough picture to extract useful information 
about the system evolution. According to the above, for the linear region we can write 
Pi{k) = A exp^r/ri) A;"' where here ni = 2, for the quintic model ti = -\/6/4, and A is a 
constant. Equating Pi{kc) = Pni{kc) at the transition we find 



kc{r) 



expi-r/nini-n)). (20) 



In practice, Eq.(20) describes the transition, defined as the intersection of the two power 
laws, remarkably well. Within the context of this model, 27r/fcc represents the size of 
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the largest clusters at that epoch. If we further assume that the intra-cluster particle 
density remains roughly constant, then we also see that the number of clusters deceases 
approximately as exp(— r/rm) where = T/(n/ — n) can be taken as a measure of the 
mean time between mergers. 

We can also obtain insight into the fractal nature of the distribution of particles in 
configuration space from this simple model. From the Weiner-Khinchine Theorem [H] , 
the correlation function ^{r) can be obtained from the Fourier transform of P(k): 

^(r) = — I dkP{k)cos{kr). (21) 



2vr _ 

Combining the contributions on either side of kc we obtain after a little algebra 

^(r) = ^Bkl^" ^J^ du u^cos{rkcu) + duu'' cos{rkcu) | . (22) 

The correlation dimension is arguably the most robust of the class of generalized fractal 
dimensions (see below). We can obtain the correlation dimension by examining the 
average distribution of particles C{R) in a ball of radius R around a given particle 



C{R) = p r dr{l + i{r)) 

J-R 



2p|i?+ -^5A;;? 



pi roc 

/ du u sin{Rkcu) + / duu^~^ sin{RkcU 



27r 

The correlation dimension is given by 

In C X 

Do = lim ; = lim ; In 

R^o In R R^o In R 




R+^ (aik'^-^'R + a^iR- 



(25) 



27r 

where, in the above, p is the average density, the second integral was expressed in terms 
of the incomplete Gamma function, only the lowest order contributions to C{R) were 
included and a; and a„/ are constants of order unity. We see that when r ^ Tm/{n + 1), 
i.e. early in the evolution, the linear term dominates and we anticipate that numerical 
analysis of the simulations will yield D2 — 1 while, for late times, the last term, 
contributed by the nonlinear regime, is dominant and D2 — —n. 

Careful observation of the clusters in the phase plane suggests a possible mechanism 
for the hierarchical, bottom-up cluster formation. The large amplitude on small scales 
initiates a large collection of small clusters. Perhaps they are seeded on the caustics 
studied by Yano and Gouda[Tl]. Once these clusters virialize, they "freeze out" of 
the expansion and become the new "particles" for the next epoch of larger clusters 
[1]. In turn, these cluster "particles" virialize and merge and the process continually 
repeats until boundary effects interfere. At large scales P{k) continues to increase as a 
result of the exponential growth in the remnants of the growing mode. On the other 
hand, at smaller scales, the freezing out of the clusters blocks the steady increase of 
P{k). Thus virialization on successively larger scales determines the shape of the power 
spectra and qualitatively explains the observed results of the simulations. Since it is 
a renormalization process, the evolution is self-similar. Of course, it is a continuous 
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process. The fact that larger "sub-clusters" become the particles for the next stage of 
virialization causes the clusters to grow in real space. In their study, Yano and Gouda 
[13] found three scaling ranges, the linear regime for large scales, a small intermediate 
range with an index of —1, and the third region which they deemed multi-caustic. 
They hypothesize that the power-law behavior is induced by the generation of caustics 
following the first crossings. The small regime, with n = —1, is associated with the 
period of single crossings. In our simulations we see hints of the intermediate regime 
in the transition regions of Fig. ([5]). In addition the index values we obtained for the 
behavior on small scales approximately correspond to theirs. Since Gabrielli et al have 
argued that the Q model and the Zeldovich approximation are equivalent , this isn't 
surprising. Additional discussion of the transition from the linear to nonlinear regime 
can be found in [ID]. 



1e+08 




Figure 5: Plots of the Power spectrum at r = 0, 2, 4, 6, 8, 10, 12, 14, 16 for the quintic 
model. Two scaling regions are observed separated by the transition from linear to 
nonlinear behavior. 
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4. Fractal Measures 

It is natural to assume that the apparently self-similar structure that develops in the 
phase plane (see figure [T]) as time evolves develops fractal geometry, but we will see that 
things are not so simple. In their earlier study of the RF model, Rouet, Feix and Jamin 
found a box counting dimension for the particle positions in /i space of about 0.75 for an 
initial water bag distribution (uniform on a rectangle in the phase plane) and a fractal 
dimension of about 0.57 in the configuration space [i.e., of the projection of the set of 
points in fi space on the position axis) [llj. As far as we know, Balian and Schaeffer were 
the first to suggest that the distribution of galaxy positions is consistent with a bifractal 
geometry [39] . Their idea was that the geometry of the galaxy distribution was different 
in the clusters and voids and, as a first approximation, this could be represented as a 
superposition of two independent fractals. Of course, their analysis was restricted solely 
to galaxy positions. Since the structures which evolve are strongly inhomogeneous, here 
we perform a multifractal analysis [43] of the configuration space. 

The multifractal formalism shares a number of features with thermodynamics 
[ISt 112]. To implement it we partitioned configuration space into cells of length /. 
At each time of observation in the simulation, a measure /ij = Ni{t)/N was assigned 
to cell i, where iVj(t) is the population of cell i at time t and is the total number of 
particles in the simulation. The generalized dimension of order q is defined by 



where Cq{l) is the effective partition function |i42j, Dq is the box counting dimension, 
Di, obtained by taking the limit g — )■ 1, is the information dimension, and D2 is the 
correlation dimension [151 As q increases above 0, the Dg provide information on 
the geometry of cells with higher population, i. e. the regions of high density or clusters. 

In practice, it is not possible to take the limit / — > with a finite sample. Instead, 
one looks for a scaling relation over a substantial range of In / with the expectation that 
a linear relation between InCg and In / occurs, suggesting power-law dependence of Cq 
on I. Then, in the most favorable case, the slope of the linear region should provide the 
correct power and, after dividing by g — 1, the generalized dimension Dg. As a rule, or 
guide, if scaling can be found either from observation or computation over three decades 
of /, then we typically infer that there is good evidence of fractal structure [B]. Also 
of interest is the global scaling index r^, where Cg ~ Z"^' for small /. From Eq. (26 we 
see that Tg and /(«) are related to each other by Dg{q — 1) = for g 7^ 1 [12]. Here we 
present the results of our fractal analysis of the particle positions on the fine (position 
only). 

If it exists, a scaling range of / is defined as the interval on which plots of In Cg 
versus In / are linear. Of course, for the special case of g = 1, we plot — S/ijln/ij vs In/ 
to obtain the information dimension [12]. If a scaling range can be found, Dg is obtained 
by taking the appropriate derivative. It is well established by proof and example that, 
for a normal, homogeneous fractal, all of the generalized dimensions are equal, while for 
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an inhomogeneous fractal, e.g., the Henon attractor, -Dq+i < Dq |13]. In the hmit of 
small /, the partition function Cg{l) can also be decomposed into a sum of contributions 
from regions of the inhomogeneous fractal sharing the same pointwise dimension a, 



where / (a) is the fractal dimension of its support [131 El 112] • Then if, for a range of q, 
a single region is dominant, we find a simple relation between the global index Tq and 

a. 



and a corresponding linear relation between In Cq and q. 

Initially the system is very cold, the fi space distribution of particles appears as 
a line, and the initial dimension is also about unity. In fact, initially the velocity of 
the particles is a perturbation. It is not large enough to allow particles to cross the 
entire system in a unit of time. After a while the fluctuations of force in the system 
destroys the approximate symmetry of the initial /i space distribution. Breaking the 
symmetry leads to the short-time dissipative mixing that results in the separation of 
the system into clusters. The behavior of the distribution of points in configuration 
space is similar. The initial dimension is nearly unity until clustering commences. At 
this time the dimensions in n space and configuration space separate. 

As time progresses, however, for the initial conditions discussed above, typically one 
dominant scaling region and a hint of a second develops. Of course this is in addition 
to the trivial scaling regions obtained for very small /, corresponding to isolated points, 
and to large / on the order of the system size, for which the matter distribution looks 
smooth. The observed size of each scaling range depended on both the elapsed time 
into the simulation and the value of q. In some instances it was possible to find good 
scaling up to 4.7 decades in / ! 

In figure joj we provide plots of In Cq versus In / in configuration space for four 
different values of q covering the range we investigated, (—4 < g < 4). To guarantee that 
the fractal structure was fully developed, we chose r = 16 for the time of observation 
and the initial conditions are those given above. There was a basic consistency in the 
results with the largest dominant scale —5.5 < InZ < 6.5 occurring for g = 4 and the 
smallest, —4 < In / < 4, for q = —4, i.e. from about 3.5 to 5 decades. In all cases there 
was a hint of a second scaling region just beyond the primary, e.g. for 6.2 < Inl < 7.6 
for the case g = 4. This scale corresponds to the turnover point in the power spectra. It 
appears to be a manifestation of the remnant of the linear region at large scales. The 
straight (dashed) lines that appear in the data are best curve fits obtained from linear 
regression for the range of the larger scaling region that is common to each of the values 
of q. Their slopes were used to compute the generalized dimensions. 

In figure [7] we illustrate the behavior of the generalized dimension Dq for the 
configuration space of the quintic model at r = 16 for —4 < g < 4. Although the 




(27) 
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Figure 6: Scaling behavior in configuration space at r = 16 for tlie quintic model. Plots 
of ^InCg versus In / are provided for four values of q: a) q=-4, b) q=0. c) q=2, d) q= 
4. The slopes have been determined using a linear regression from points 12 to 20. 



embedding dimension is c? = 1, Dq is about 0.65 so the distribution is definitely fractal. 
For g < , the curve is an increasing function. This unphysical behavior suggests that 
there is simply a lack of sufficient data in the low-density regions to provide reliable 
dimensions. On the other hand, for g > , we see the decreasing behavior characteristic 
of a multifractal. In figure [8] we present a plot of vs q constructed from the same 
data. Note the constant value for g < 0. Comparison with equation (28) suggests that 
the pointwise dimension vanishes in the under-dense regions and that there is simply 
insufficient data to fix the type of geometry [12]. This is consistent with figure [7j For 
the case where q > 0, the nearly linear behavior suggests that equation (28) applies here 
as well and yields a pointwise dimension a of about 0.65, and a support with Hausdorff 
dimension f{a) of about 0.7. It is not immediately obvious how to reconcile figures [T] 
and [s] for g > but it is interesting that both a and f{a) are very close to the value 
of Dq. In varying the initial conditions by selecting different random number seeds, 
there are slight variations in the overall picture. For g > we find Dg consistently 
in the neighborhood [0.62, 0.65] but the degree of the multi-fractal nature, gauged by 
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the negative slope in figure [7| varies somewhat. Perhaps the multifractal character 
suggested by the decrease of D{q) with increasing q is an artefact due to the increasing 
sparseness of denser and denser regions. However, since simulations of the RF model 
exhibit similar structure formation and a flat D{q), this is doubtful. 




Figure 7: Generalized dimension Dg vs q in configuration space for the quintic model 
at r = 16. 

Part of our goal is to compare how fractal geometry arises in a family of related 
models. So far we have presented results for the quintic, or Q, model. However we 
have also carried out similar studies of the Rouet-Feix (RF) model and the Hamiltonian 
model without friction (which can be obtained from either of the former by nullifying 
the first derivative contribution in equation (16)). In table [l] we compare the correlation 



dimension D2 and the exponents generated by the power spectrum and correlation 
function for all three models at scaled times r = 13 and 16, where the clustering is 
highly developed. Writing P{k) ~ and ^(r) ~ r~'^, in principle there are simple 
relations between D2, n, and 7 ; namely D2 = 1 — '~f = —n. The first equality can 
be easily derived from the definition of correlation dimension (see e.g. [HUl |12]and the 
discussion in the previous section), while the second is due to the Wiener- Khintchine 
theorem, i.e. the power spectrum is the Fourier transform of the correlation function 
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-0.5 



Figure 8: Global scaling index Tq vs. q in configuration space for the quintic model at 
r = 16. 



[S]. We see from table [Tjthat these relations are approximately obeyed in these systems 
with 2^^ particles. The agreement is superior for the Q and RF models. Experience 
has shown that, for them, as the system size is increased, so also is the agreement. 
While there are similarities in the fractal structure of the RF and Q models, they are 
not the same. For example in the quintic model the generalized fractal dimension Dg is 
consistently greater than the corresponding dimension in the RF model showing that the 
inhomogeneity is stronger in the latter although the larger friction constant of the RF 
model causes clustering to commence later in the system evolution. Another difference 
is that plots of Dqvsq for positive q for the RF model are virtually constant, suggesting 
monofractal behavior. The values obtained for n are consistent with those found in the 
simulations carried out by Yano and Gouda using the Zeldovich approximation in the 
"multicaustic" regime [TB] . 

The Hamiltonian system was also investigated by Gabrielli et al |23] . Of the three 
models, it exhibits the most rapid cluster formation. However, in other aspects the 
behavior differs sharply from the dissipative systems. First, the scaling behavior is less 
robust. As shown in figures |9| plots of the power spectra at T = 13, 16 have a smaller 
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scaling range. In addition, the fractal properties are much closer to an ordinary system. 
In figure 10 we see that the generalized dimensions Dg 0.9, much greater than the RF 
model, which presents the greatest structure with Dq ~ 0.5. As can be seen in table 
[1} the connection between the exponents is weaker, no doubt due to the smaller scaling 
range. The relation between the exponents is only exact if the correlation function 
and power spectra are power-laws. Finally, if we extend the simulations beyond the 
time where boundary conditions can be neglected, say for T > 16, we see a remarkable 
difference. For such times the matter distribution in the dissipative systems eventually 



collapses to just a single cluster. In contrast, as we see in figure 11, the Hamiltonian 
system becomes space filling. It appears to revert to the stationary state for this system 
introduced by Valageas [HI [19] . 



p(k) 

P(k) averaged on 41 points 
slope n=-1 ,05 




P(k) 

P(k) averaged on 41 points 
slope n=-.88 




(b) 

Figure 9: The Power spectrum for the Hamiltonian model : (a) at r = 13 it exhibits 
a scaling region in k from about 0.01 to 1, with a scaling exponent —n = 1.05 (b) at 
r = 16, the scaling range is [0.01;2], and —n = 0.88. 



There is a sensitive dependence of the evolution on the initial conditions and, 
in particular, on the initial spectral index n. Although the cluster distributions in 
configuration space are similar, more information is available in the phase plane. We 
find that the smaller the initial index, the greater is the predominance of large clusters 



in the initial state, and these are very unstable. In figures 12 and 13 we show snapshots 
of the evolution of the Q model with n = 1 and 0. We see from the figures that there 
is a strong correlation between the initial condition and the type of structure formation 
that develops. For the case where the initial particle distribution is white noise (n = 0), 
examine the /i space snapshot at the time r = 13. Note that the system has not broken 
apart as we would expect from the simulations with n = 2. Moreover, the characteristic 
bottom-up behavior, where the smaller clusters combine to form larger ones as time 
evolves, is missing. Instead there is a single, large, connected, structure with a few 
small "subclusters" distributed within the pattern. Clearly hierarchical clustering is 
absent. Examining the snapshots for n = 1 we see that the behavior is intermediate. 
Both small and large clusters are present, but the bottom-up scenario is not nearly as 



Cosmology in One Dimension: Fractal Geometry, Power Spectra and Correlation 25 



0.9 
0.8 
0.7 
0.6 

cr 
Q 

0.5 
0.4 
0.3 
0.2 

-4-3-2-1 1 2 3 4 

q 



Figure 10: Generalized dimension Dq vs q in configuration space for tlie Hamiltonian 
model at r = 16. 

robust as for n = 2. A clue for the origin of this behavior may be found in the work 
of Gabrielli et al mentioned earlier |2H [25]. In their investigation of the statistical 
properties of static distributions of infinite systems that are on average uniform, they 
found that the variance of the force diverges if the structure factor is characterized by 
white noise {n = 0) and is marginal for n = 1. Alternatively, the source of the behavior 
may arise from dynamical considerations. Energy at small scales is required to establish 
the presence of small clusters that can initiate a process of hierarchical virialization. 
For n = there is simply too much competition from large clusters and the hierarchical 
process is broken. Although the hierarchical process appears to start in a few locations, 
it doesn't have sufficient time to develop. Long wavelengths drive the overall simulation. 

5. Summary and Conclusions 

Here and in our earlier work we have seen that hierarchical clustering is an intricate, 
highly choreographed dance in phase space that requires an initially cold and uniform 
system. The three models that were studied here differ only in the value of the friction 
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Figure 11: Evolution in configuration and (i space for the Hamiltonian model with 2^^ 
particles from r = to r = 20. The initial distribution is such that the density power 
spectrum has a power-law of index n = 2. 

constant, a: a =^)^) , respectively, for the RF, Q, and Hamiltonian models. The Q 
(Zeldovich) model has a hybrid geometry that may actually represent the local behavior 
of a 3+1 dimensional system for a short time, but cannot reproduce its real clustering 
properties in the highly nonlinear regime. Alternatively, the RF model is completely 
self-consistent with the transformation to comoving coordinates in 1+1 dimensions. It 
is useful for unambiguously studying nonlinear behavior, but has a weaker connection to 
the 3+1 dimensional universe. In fact, simulations of both the Q and RF models show 
remarkably similar behavior. The Hamiltonian model provides a contrasting dynamical 
system with no dissipation in the comoving coordinate system. Like the Q and RF 
models, it undergoes rapid hierarchical clustering for a short time span. However, unlike 
them, on long time scales the system becomes space filling and appears to revert to the 
stable stationary state investigated by Valageas [iHl [19] . 

To mimic cosmological investigations, the initial positions were selected by sampling 
a power-law spectrum with random phases. The velocities were chosen to insure that 
the system is in the growing mode in phase space. The three models were studied 
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Figure 12: Evolution in configuration and fi space for tlie quintic model witli 2^^ 
particles from r = to r = 16. The initial distribution is such that the density power 
spectrum has a power-law of index n = 1. 
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Figure 13: Evolution in configuration and fi space for the quintic model with 2^^ 
particles from r = to r = 13. The initial distribution is such that the density power 
spectrum has a power-law of index n = 0. 



with a range of initial power-law indices n. Qualitative behavior was observed for 
n = 3,4 while extensive studies were carried out for n = 0, 1, 2. Both the RF and 
Q models exhibited very similar behavior. For the case n = 2, initially the systems 
remain fairly uniform. After some time, around r = 8 for the Q model, the systems 
break into small clusters. As time progresses further the clusters coalesce into larger 
ones, supporting the well known bottom-up scenario of hierarchical cluster formation 



Cosmology in One Dimension: Fractal Geometry, Power Spectra and Correlation 28 

predicted for the real universe. Concurrently the variance of the velocities grows 
approximately exponentially. The power spectrum was monitored at various times and 
continued to show scale-free, power-law, behavior. In the strongly clustered epoch it 
turns around at small scales to a negative value supporting the observed structure of 
clusters of increasing size. It was shown that a simple analytic representation of the 
power spectra captured the main features of the evolution, including the correct time 
dependence of the crossover from the linear to nonlinear regime and the transition from 
regular to fractal geometry. The underlying mechanism for the observed behavior may 
be a process of "hierarchical virialization" whereby virialized clusters that have "frozen 
out" from the cosmic expansion become the new "particles" for the next virialization 
sequence. The measurements were consistent with the self-similar evolution found in 
the simulations of Yano and Gouda using the Zeldovich model |13j . 

In the well clustered regime, a multifractal analysis of the configuration space was 
performed. For the RF and Q models, in contrast with our earlier work with isothermal 
and waterbag initial conditions [HI [17], here a single dominant scaling range of the 
partition function Cq was found for all values of q considered. Generalized dimensions 
Dg and scaling index Tg were evaluated for —4 < q < 4. Robust fractal behavior 
was inferred for g > from the values of Dq. For the Q model the gradual decrease 
in Dq for g > suggests that the system is also weakly multifractal. The RF model 
exhibited the greatest inhomogeneity. It appeared nearly monofractal with D2 = 0.5. In 
addition, the two-body correlation function was constructed and also exhibited a power- 
law dependence on the displacement. For large samples, on the order of 10^ or more 
particles, the mathematical relation between the correlation dimension D2, the spectral 
index n and the correlation function decay exponent 7, D2 = —n = 1 — 7, was nicely 
satisfied. For the Q model at r = 13, D2 = 0.63 over a finite scale range. For n > 2 
there was evidence of hierarchical clustering whereas for n < 2 there was a tendency to 
form a few large clusters, but in the absence of self-similar, bottom-up, evolution. Like 
the Q and RF models, the Hamiltonian version undergoes rapid hierarchical clustering 
for a short time span. However, the scaling behavior was less robust, and the agreement 
among the exponents was weaker. In conclusion, one-dimensional models can exhibit 
the type of evolution characteristic of the real universe. Their evolution does reveal 
fractal behavior, but on a limited length scale that changes with time. It also shows 
that hierarchical structure formation depends sensitively on the initial conditions. Their 
three dimensional counterparts are useful for constraining models of the early universe. 

A number of interesting questions about these models still remain. In forthcoming 
work we will explore the connection between symmetry and the various one-dimensional 
models and show how additional models analogous to the Q type can be constructed. We 
will also establish a connection between the increase in the variance of the velocity and 
virialization. While the fractal properties of the over-dense clusters appear to be well 
represented in the current work, the geometry of the under-dense "voids" has not been 
revealed. This is under investigation with a greater number of particles per cluster to 
improve the resolution. A complementary approach that may retain more information in 
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the under-dense regions is provided by the numerical integration of the Vlasov equation. 
To finally address the conjecture of Balian and Schaeffer concerning the bifractality of 
the universe for this class of models, this must be completed. The approximate analytic 
representation of the power spectra introduced in Section (III) suggests that the low 
density inter-cluster regions have fractal dimension unity, but these results must be 
considered preliminary. A compelling issue is the connection between cosmologies in 
different dimensions. In an earlier work we showed that the RF model was able to 
crudely reproduce the time of galaxy formation in a matter-dominated universe [I6] . It 
would be intriguing to be able to establish a more precise correspondence. 
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